imports
In [1]:
import plotly.offline as py
import plotly.graph_objs as go
py.init_notebook_mode(connected=True)
#numpy for calculations
import numpy as np
In [2]:
"""utils.py"""
import numpy as np
def mesh2d(xlim, ylim, n=5):
"""Create 2d mesh in sepecifies x and y axes limits, with number of points for every dimension separate."""
if isinstance(n, int):
xx = np.linspace(xlim[0],xlim[1],n)
yy = np.linspace(ylim[0],ylim[1],n)
elif isinstance(n, list):
xx = np.linspace(xlim[0],xlim[1],n[0])
yy = np.linspace(ylim[0],ylim[1],n[1])
else:
raise Exception("Wrong number of points parameter")
return np.meshgrid(xx, yy)
In [3]:
"""utils.py"""
def normalize(v):
"""Normalizes a 3d vector v, returns a 3d vector."""
magnitude = np.sqrt(v[0]**2+v[1]**2+v[2]**2)
if magnitude==0:
raise ValueError("Zero vector cannot be normalized.")
else:
return v/magnitude
In [4]:
"""utils.py"""
def jsonify(obj):
"""Transforms a graphic object (go) into json data dictionary.
@author Nick Metelski
@since 27.07.17"""
go_obj = obj.goify()
json_obj = dict()
json_obj['name'] = go_obj['name']
json_obj['x'] = list(go_obj['x'])
json_obj['y'] = list(go_obj['y'])
json_obj['z'] = list(go_obj['z'])
return json_obj
We need classes for Point, Line and Plane for uniform output:
In [5]:
"""point.py"""
class Point:
"""Point class to make returns from intersections more reasonable.
@author Nick
@since 25.07.17"""
def __init__(self,position):
self.pos = np.array(position)
def getXYZ(self, layout=None):
"""Generate x,y,z data based on the layout box. Returns tuple (x,y,z).
@author Nick Metelski
@since 27.07.17"""
return np.array([self.pos[0]]), np.array([self.pos[1]]), np.array([self.pos[2]])
def goify(self):
"""Transform a point into graphics object (go).
@author Nick Metelski
@since 25.07.17"""
pt = go.Scatter3d(
mode="markers",
x=[self.pos[0]],
y=[self.pos[1]],
z=[self.pos[2]]
)
return pt
Testing
In [6]:
pt = Point([1,2,3])
print(pt.goify())
print(jsonify(pt))
In [7]:
"""line.py"""
import plotly.graph_objs as go
import numpy as np
class Line:
"""Line class for intersections simulation
@author Nick Metelski
@since 25.07.17"""
def __init__(self, vec, offset):
self.vec = normalize(np.array(vec)) #normalize direction vector
self.offset = np.array(offset) #cast to numpy array
#minimize the offset to be minimum distance
self.offset = self.offset - (np.dot(self.offset,self.vec) * self.vec)
def getXYZ(self, layout=None):
"""Generate x,y,z data based on the layout box. Returns tuple (x,y,z).
@author Nick Metelski
@since 26.07.17"""
#TODO import data from layout
t = np.linspace(0,1,2) #the parameter
xx = self.offset[0]+t*self.vec[0] #generate xpoints
yy = self.offset[1]+t*self.vec[1] #for y points
zz = self.offset[2]+t*self.vec[2] #for z points
return xx, yy, zz
def goify(self,layout=None):
"""Export the line into graphics object
@author Nick Metelski
@since 26.07.17"""
xx, yy, zz = self.getXYZ()
line = go.Scatter3d(
mode="lines",
x=list(xx),
y=list(yy),
z=list(zz),
line = dict(
color = ('rgb(205, 12, 24)'),
width = 10)
)
return line
Testing
In [8]:
line = Line([1,2,3],[1,1,1])
print(line.goify())
print(jsonify(line))
Equation of a plane: $$\vec{n} \cdot (\vec{r} - \vec{r_0}) = 0$$
In [9]:
"""plane.py"""
"""Plane class for intersections simulation
@author Nick
@since 25.07.17"""
import plotly.graph_objs as go
#import .line #import local line module
import numpy as np
class Plane:
"""Planes are defined by their normal and offset vector."""
def __init__(self, normal, offset):
self.normal = normalize(np.array(normal)) #normalize normal vector
self.offset = np.array(offset)
self.offset = np.dot(self.offset,self.normal)*self.normal
def getXYZ(self, xlim=[-1,1], ylim=[-1,1], zlim=[-1,1], n=2):
"""Generate x,y,z data based on the layout box. Returns tuple (x,y,z).
@author Nick Metelski
@since 25.07.17"""
if self.normal[2] == 0: #check if z is zero, then we have to generate x or y from other meshes
if self.normal[1] == 0: #check if z and y is zero, then we have to generate x from yz mesh
if self.normal[0] == 0:
return ValueError("Normal vector is zero vector.")
else:
#cannot generate z but can y, try generating y for xz mesh
y, z = mesh2d(ylim, zlim, n)
x = (np.dot(self.normal,self.offset)-self.normal[1]*y-self.normal[2]*z)/self.normal[0]
else:
#cannot generate z but can y, try generating y for xz mesh
#self.normal[2] = 0.01 # TODO THIS IS VERY CRUDE
x, z = mesh2d(xlim, zlim, n)
y = ((np.dot(self.normal,self.offset)
- self.normal[0]*x
- self.normal[2]*z)
/ self.normal[1])
else:
#try generating z
x, y = mesh2d(xlim, ylim, n)
#Generate plane z-values array
z = (np.dot(self.normal,self.offset)-self.normal[0]*x-self.normal[1]*y)/self.normal[2]
return [list(i) for i in x],[list(i) for i in y],[list(i) for i in z]
def goify(self, layout=None):
"""Export the plane into graphics object.
@author Nick Metelski
@since 25.07.17"""
xx,yy,zz = self.getXYZ()
surf = go.Surface(
x=xx,
y=yy,
z=zz
)
return surf
Testing
In [10]:
plane = Plane([1,2,3],[1,1,1])
print(plane.goify())
Generate some objects
In [11]:
point1 = Point([0.5,0.0,0.2])
point2 = Point([0.1,-0.4,0.3])
line1 = Line([2,1,0],[0,0.1,0])
line2 = Line([2,1,1],[0.2,0.1,0])
plane1 = Plane([1,1,1],[0,0.1,0])
plane2 = Plane([2,0,-1],[-0.2,0.1,0.3])
objlist = [point1, point2, line1, line2, plane1, plane2]
data = [obj.goify() for obj in objlist]
Define layout
In [12]:
u
Plot the figure
In [13]:
fig=go.Figure(data=data,layout=layout)
py.iplot(fig)
We need points as well as planes and lines.
In [14]:
"""algebra.py"""
"""This contains algebra methods to calculate intersections etc."""
def intersection(obj1, obj2):
"""Intersection between two algebra 3d objects.
@author Nick Metelski
@since 26.07.17"""
#plane-plane
if isinstance(obj1,Plane) and isinstance(obj2,Plane):
return _plane_plane_intersection(obj1,obj2)
#plane-line
elif isinstance(obj1,Plane) and isinstance(obj2,Line):
return _plane_line_intersection(obj1,obj2)
elif isinstance(obj1,Line) and isinstance(obj2,Plane):
return _plane_line_intersection(obj2,obj1)
#plane-point
elif isinstance(obj1,Plane) and isinstance(obj2,Point):
return _plane_point_intersection(obj1,obj2)
elif isinstance(obj1,Point) and isinstance(obj2,Plane):
return _plane_point_intersection(obj2,obj1)
#line-line
elif isinstance(obj1,Line) and isinstance(obj2,Line):
return _line_line_intersection(obj1,obj2)
#line-point
elif isinstance(obj1,Line) and isinstance(obj2,Point):
return _line_point_intersection(obj1,obj2)
elif isinstance(obj1,Point) and isinstance(obj2,Line):
return _line_point_intersection(obj2,obj1)
#point-point
elif isinstance(obj1,Point) and isinstance(obj2,Point):
return _point_point_intersection(obj1,obj2)
#wrong params
else:
raise TypeError("Invalid parameter types - please pass intersections.Point, intersections.Line or intersections.Plane.")
def _plane_plane_intersection(p1,p2):
"""Private function; do not use. Use intersection(obj1,obj2) instead.
plane-plane intersection submethod. Raises exceptions or returns Line.
@author Nick Metelski
@since 26.07.17"""
#cross-product
cross = np.cross(p1.normal,p2.normal)
if np.all(cross==0):
#planes are parallel or overlap
raise ArithmeticError("Planes are parallel.")
else:
#sample point: x=0
mat = [[p1.normal[1],p1.normal[2]],[p2.normal[1], p2.normal[2]]]
axis = 0 #keep track of which axis we zero out; 0=x, 1=y, 2=z
#NOTE: linalg.solve can raise np.linalg.LinAlgError exception when x=0 results in singular matrix
#we have to check some other cases
#premise: the line has to intersect at least one of the planes: xy, yz or xz.
if np.linalg.matrix_rank(mat) == 1:
mat = [[p1.normal[0],p1.normal[2]],[p2.normal[0], p2.normal[2]]]
axis = 1
if np.linalg.matrix_rank(mat) == 1:
mat = [[p1.normal[0],p1.normal[1]],[p1.normal[0], p2.normal[1]]]
axis = 2
rhs = [np.dot(p1.normal,p1.offset),np.dot(p2.normal,p2.offset)]
sol = np.linalg.solve(mat,rhs)
if axis == 0:
return Line(cross, [0,sol[0],sol[1]])
if axis == 1:
return Line(cross, [sol[0],0,sol[1]])
if axis == 2:
return Line(cross, [sol[0],sol[1],0])
def _plane_line_intersection(plane,line):
"""Private function; do not use. Use intersection(obj1,obj2) instead.
Intersection of a line an a plane. Raises exceptions or returns Point.
@author Nick Metelski
@since 26.07.17"""
check = np.dot(plane.normal,line.vec)
if np.all(check==0):
#plane and line are parallel or overlap
raise ArithmeticError("Plane and line are parallel.")
else:
#there is an explicit formula for parameter of the line for which intersection is met:
#$$t = \frac{\vec{n} \cdot (\vec{d_p}-\vec{d_v})}{\vec{n} \cdot \vec{v}}$$,
#where n is normal to plane, v is line's direction, rp is plane offset, rv is vector offset
t = np.dot(plane.normal,(plane.offset-line.offset))/np.dot(plane.normal,line.vec)
return Point(line.offset+t*line.vec)
def _plane_point_intersection(plane,point):
raise NotImplementedError("Plane-point intersections not implemented.")
def _line_line_intersection(line1,line2):
raise NotImplementedError("Line-line intersections not implemented.")
def _line_point_intersection(line,point):
"""Private function; do not use.
Intersection of a line and point. Intersection is a Point.
Raises exceptions or returns Point.
@author Nick Metelski
@since 27.07.17"""
offset_diff = point.pos - line.offset #difference in offsets must be parallel to direction vector
cross = np.cross(line.vec,offset_diff) # this is for checking if there exists an intersection
if cross == 0:
return Point(point.pos) #return a point coinciding with the original point passed as param
else:
raise ArithmeticError("Point does not lie on the line.")
def _point_point_intersection(point1, point2):
"""Private function; do not use. Use intersection(obj1,obj2) instead.
Intersection of a two points. Intersection is the point if the point are the same position.
Raises exceptions or returns Point.
@author Nick Metelski
@since 27.07.17"""
if np.all_equal(point1,point2):
#just return one of the moints when they coincide
return Point(point1.pos)
else:
raise ArithmeticError("Points do not coincide.")
In [15]:
plane1 = Plane([1,1,1],[0,0,0])
plane2 = Plane([1,-1,1],[0,0,0])
plane3 = Plane([-1,1,0],[0,0.5,0])
line1 = intersection(plane1,plane2)
line2 = intersection(plane2,plane3)
objs = [plane1, plane2, plane3, line1, line2]
data = [obj.goify() for obj in objs]
fig=go.Figure(data=data,layout=layout)
py.iplot(fig)
In [16]:
try:
plane = Plane([1,0,0],[0,0.5,0])
line = Line([0,1,1],[-.5,0,0.5])
intersec = intersection(plane,line)
fig=go.Figure(data=[plane.goify(),line.goify(),intersec.goify()],layout=layout)
py.iplot(fig)
except:
raise
In [ ]:
help(Point)
help(_point_point_intersection)
In [ ]: